#Base plate pull-up model based on Malhotra and Veletsos (1994) as implemented in Bakalis et al. (2017)
 
import openseespy.opensees as op
import opsvis as opv
import os
import numpy as np
import matplotlib.pyplot as plt
import time
import Units as ut

op.wipe();

print('Model generation started...')

direc = './ResultsBP'

if(os.path.exists(direc)=='False'):
    os.mkdir(direc);


R = 23.47*ut.m #Tank Radius in m
H = 19.95*ut.m #Contained Liquid Height in m
n = 8 #Number of beams used to modelling the base plate, Bakalis et al. (2017)
tw = 22.23*ut.mm #wall thickness in mm
tb = 6.4*ut.mm #base plate thickness
ta = 10*ut.mm #annular ring thickness
ro = 998*ut.kg_m3 #liquid density in kg/m3  
Ew = 1000*ut.MPa #Young Modulus of rigid foundation in MPa 
Es = 210000*ut.MPa #Young Modulus of steel in MPa
Fys = 337*ut.MPa #Yield strength of steel, MPa
epy = Fys/Es
Bs = 0.01 #strain hardening ratio
v = 0.3 #Poisson ratio of steel

fillr = 0.95 #Pencentage of tank's height filled with liquid 

bw = 2*np.pi*R/n  #uniform width, Bakalis et al. (2017)
ros = 7850*ut.kg_m3 #density of steel
Ww = ros*ut.g*H*bw*tw #weight of wall over one side of the strip
mr = 35*ut.ton
Wr = mr*ut.g #weight of roof
Wrr = Wr/n
ph = ro*ut.g*H*fillr  #Hydrostatic pressure at tank bottom

ktt = Es*bw*tw**2*(tw/R)**0.5/(2*(3*(1-v**2))**0.75) #Rotational stiffness due to wall, Bakalis et al. (2017)
kuu = Es*bw*(tw/R)**1.5/((3*(1-v**2))**0.25)  #Axial stiffness due to wall, Bakalis et al. (2017) 
Mr = bw*R*tw*ph/(2*(3*(1-v**2))**0.5)   #Point moment at bas plate wall junction due to hydrostatic pressure, Bakalis et al. (2017)
Nr = bw*(R*tw)**0.5*ph/(2*(3*(1-v**2))**0.25) #Point radial load due to hydrostatic pressure, Bakalis et al. (2017) 

kw = Ew*bw #stiffness of zero-length Winkler spring at tank foundation interface

lb = 50*ut.mm  #length of beam sub-element aprox. 15tb, Bakalis et al. (2017) 
nb = int(2*R/lb)   #number of beam sub-elements

la = 500*ut.mm    #length of annular ring 
na = int(round(la/lb,0))   #number of elements in the annular rings 
wph = ph*bw  #distributed load on beam due to hydrostatic pressure
print(wph)
#Total Weight of Liquid
Wtot = wph*2*R
# Point Load in Nodes
wpn = Wtot/(nb+1)


op.model('basic', '-ndm', 2, '-ndf', 3) 

#Define nodes
nodeTag1 = 10000
nodeTag2 = 20000

for i in range(nb+1):
    op.node(nodeTag1+i,i*lb,0) 
 

#op.fix(nodeTag1+nb,1,0,1)    

for i in range(1,nb):
    op.node(nodeTag2+i,i*lb,0)
    op.fix(nodeTag2+i,1,1,1)
    
    
op.fix(nodeTag1,1,1,0) #pined support at compression edge

nodeMaster = 30000

op.node(nodeMaster,nb*lb,0)
op.fix(nodeMaster,1,0,1)

#Define Materials
steelMatTag = 1
op.uniaxialMaterial('Steel01', steelMatTag, Fys, Es, Bs)
winkMatTag = 2
op.uniaxialMaterial('ENT',winkMatTag,kw)
kttMatTag = 3
op.uniaxialMaterial('Elastic',kttMatTag,ktt)
kuuMatTag = 4
op.uniaxialMaterial('Elastic',kuuMatTag,kuu)
rigidMatTag = 5
op.uniaxialMaterial('Elastic',rigidMatTag,ut.Ubig)
flexMatTag = 6
op.uniaxialMaterial('Elastic',flexMatTag,ut.Usmall)

#Create base plate fiber section
secTag = 1
op.section('Fiber',secTag)
# number of fibers in z-direction
nfZ = 4
# number of fibers in y-direction
nfY = int(round(bw/5,0))
y1 = -bw/2
yJ = bw/2
z1 = -tb/2
zJ = tb/2
op.patch('rect', steelMatTag, nfZ, nfY, z1, y1, zJ,yJ)
#Create annular ring fiber section
secTaga = 2
op.section('Fiber',secTaga)
# number of fibers in z-direction
y1a = -bw/2
yJa = bw/2
z1a = -ta/2
zJa = ta/2
op.patch('rect', steelMatTag, nfZ, nfY, z1a, y1a, zJa,yJa)

#Create elements

#Winkler Springs
eleTagW = 20000
for i in range(1,nb):
    op.element('zeroLength',eleTagW+(i-1),nodeTag2+i,nodeTag1+i,'-mat',winkMatTag,'-dir',2)

#Beam Elements
intgTag = 1
intgTaga = 2
numInt = 10
maxIter = 10
tol = 1e-12
#Define geometric transformation
BPTransf = 'Corotational'
transfTag = 1
op.geomTransf(BPTransf,transfTag)

op.beamIntegration('Legendre',intgTag,secTag,numInt)
op.beamIntegration('Legendre',intgTaga,secTaga,numInt)
eleTagBeam = 30000
for i in range(0,nb):
    if (i<na or i>=nb-na):
        op.element('forceBeamColumn',eleTagBeam+i,nodeTag1+i,nodeTag1+(i+1),transfTag,intgTaga,'-iter',maxIter,tol)
    else:
        op.element('forceBeamColumn',eleTagBeam+i,nodeTag1+i,nodeTag1+(i+1),transfTag,intgTag,'-iter',maxIter,tol)
ktueleTag = 1
op.element('zeroLength',ktueleTag,nodeMaster,nodeTag1+nb,'-mat',kuuMatTag,rigidMatTag,kttMatTag,'-dir',1,2,3)


#Run Gravity Analysis
constTS = 1
op.timeSeries('Constant', constTS)
patternTagG = 1
op.pattern('Plain',patternTagG,constTS)

for i in range(nb+1):
      op.load(nodeTag1+i,0,-wpn,0)
#op.eleLoad('-range',eleTagBeam,eleTagBeam+nb-1,'-type','-beamUniform',-wph)
op.load(nodeTag1+nb,Nr,0,-Mr)
#op.load(nodeTag1+nb,0,-(Ww+Wrr),0)

Tol = 1e-8 # convergence tolerance for test
NstepGravity = 10
DGravity = 1/NstepGravity
op.integrator('LoadControl', DGravity) # determine the next time step for an analysis
op.numberer('RCM') # renumber dof's to minimize band-width (optimization), if you want to
op.system('SparseGeneral', '-piv') # how to store and solve the system of equations in the analysis
op.constraints('Plain') # how it handles boundary conditions
op.test('EnergyIncr', Tol, 50) # determine if convergence has been achieved at the end of an iteration step
op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
op.analysis('Static') # define type of analysis static or transient
op.analyze(NstepGravity) # apply gravity
op.loadConst('-time', 0.0) #maintain constant gravity loads and reset time to zero
print('Static analysis completed')



#Define recorders
op.recorder('Node', '-file', direc+'/disp_Cyclic.out','-time', '-node', nodeMaster, '-dof', 2, 'disp')


#Run Cyclic Pushover Analysis
print('Running Cyclic Pushover...')
#displist = [50*ut.mm,0,100*ut.mm,0,150*ut.mm,0,200*ut.mm,0,250*ut.mm,0,300*ut.mm,0,500*ut.mm,0,1000*ut.mm,0]
displist = [3*ut.mm,0,3*ut.mm,0,10*ut.mm,0,10*ut.mm,0,20*ut.mm,0,20*ut.mm,0,30*ut.mm,0,30*ut.mm,0,50*ut.mm,0,50*ut.mm,0,100*ut.mm,0,100*ut.mm,0,150*ut.mm,0,150*ut.mm,0,200*ut.mm,0,200*ut.mm,0,300*ut.mm,0,300*ut.mm,0,500*ut.mm,0,500*ut.mm,1000*ut.mm,0,1000*ut.mm,0]
dstep = 0.1*ut.mm
IDctrlNode = nodeMaster			# node where disp is read for disp control
IDctrlDOF  = 2			        # degree of freedom read for disp control (1 = x displacement)

testParams = [1.e-4, 50]

op.numberer('RCM') # renumber dof's to minimize band-width (optimization), if you want to
op.system('BandGeneral') # how to store and solve the system of equations in the analysis
op.constraints('Plain') # how it handles boundary conditions
op.test('NormDispIncr',*testParams) # determine if convergence has been achieved at the end of an iteration step
op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
op.analysis('Static') # define type of analysis static or transient

POtag = 100;
linTS = 2
op.timeSeries('Linear',linTS)

op.pattern('Plain', POtag, linTS)     
op.load(IDctrlNode,0,1,0)

ok = 0
for j in range(0, len(displist)):
	if ok != 0:
		break
	dispnew = displist[j]
	if j > 0:
		dispold = displist[j-1]
	else:
		dispold = 0
	disp = dispnew - dispold
	if disp > 0:
		dsteps = dstep
	else:
		dsteps = -dstep
	nstep = int(disp/dsteps)
	op.integrator('DisplacementControl', IDctrlNode, IDctrlDOF, dsteps) # determine the next time step for an analysis
	for i in range(nstep):
		ok = op.analyze(1)
		
		if ok != 0:
			testParams2 = [1.e-4, 2000]
			op.test('EnergyIncr',*testParams2)
			op.algorithm('Newton','-initial')
			print("Trying Newton with Initial Tangent ..")
			ok = op.analyze(1)
			op.test('NormDispIncr',*testParams) 
			op.algorithm('Newton') 
		
		if ok != 0:
			op.algorithm('Broyden',50)
			op.test('RelativeNormDispIncr',*testParams2)
			print("Trying Broyden ..")
			ok = op.analyze(1) 
			op.test('NormDispIncr',*testParams) 
			op.algorithm('Newton') 
			
		if ok != 0:
			op.algorithm('NewtonLineSearch')
			print("Trying NewtonWithLineSearch ..")
			ok = op.analyze(1)
			op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
			
		if ok == 0:
			dsp = round(op.nodeDisp(IDctrlNode,IDctrlDOF),3)
			print(f'Displacement {dsp} reached of {dispnew}')
		else:
			break

if(ok != 0):
	print('Problem with Pushover')
else:
	print('Done')
		




